function stability_analysis()
% STABILITY_ANALYSIS LU分解数值稳定性分析
% 
% 分析LU分解在不同条件下的数值稳定性

fprintf('=== LU分解数值稳定性分析 ===\n\n');

%% 1. 选主元的重要性
fprintf('1. 选主元的重要性\n');
fprintf('目标：比较有无选主元的LU分解稳定性\n\n');

% 构造需要选主元的矩阵
A1 = [1e-20, 1; 1, 1];
fprintf('测试矩阵 A1 =\n');
disp(A1);
fprintf('条件数: %.2e\n', cond(A1));

% 不选主元的LU分解
try
    [L1_no_pivot, U1_no_pivot] = lu_no_pivot(A1);
    fprintf('\n不选主元LU分解:\n');
    fprintf('L1 =\n'); disp(L1_no_pivot);
    fprintf('U1 =\n'); disp(U1_no_pivot);
    
    % 验证分解
    error1_no_pivot = norm(A1 - L1_no_pivot * U1_no_pivot, 'fro');
    fprintf('重构误差: %.2e\n', error1_no_pivot);
    
    % 增长因子
    growth_factor_no_pivot = max(abs(U1_no_pivot(:))) / max(abs(A1(:)));
    fprintf('增长因子: %.2e\n', growth_factor_no_pivot);
    
catch ME
    fprintf('不选主元LU分解失败: %s\n', ME.message);
end

% 选主元的LU分解
[L1_pivot, U1_pivot, P1] = lu(A1);
fprintf('\n选主元LU分解:\n');
fprintf('P1 =\n'); disp(P1);
fprintf('L1 =\n'); disp(L1_pivot);
fprintf('U1 =\n'); disp(U1_pivot);

% 验证分解
error1_pivot = norm(P1*A1 - L1_pivot * U1_pivot, 'fro');
fprintf('重构误差: %.2e\n', error1_pivot);

% 增长因子
growth_factor_pivot = max(abs(U1_pivot(:))) / max(abs(A1(:)));
fprintf('增长因子: %.2e\n', growth_factor_pivot);

%% 2. 增长因子分析
fprintf('\n2. 增长因子分析\n');
fprintf('目标：分析不同矩阵的增长因子\n\n');

% 测试不同类型的矩阵
test_matrices = {
    randn(10, 10),                    % 随机矩阵
    hilb(8),                          % Hilbert矩阵
    gallery('clement', 10),           % Clement矩阵
    gallery('kms', 10, 0.5),         % KMS矩阵
    triu(randn(10, 10)) + eye(10)     % 上三角矩阵
};

matrix_names = {'随机矩阵', 'Hilbert矩阵', 'Clement矩阵', 'KMS矩阵', '上三角矩阵'};

fprintf('增长因子比较:\n');
fprintf('矩阵类型        尺寸    条件数        无主元增长因子    有主元增长因子\n');
fprintf('------------    ----    ------        --------------    --------------\n');

for i = 1:length(test_matrices)
    A = test_matrices{i};
    n = size(A, 1);
    cond_A = cond(A);
    
    % 不选主元（如果可能）
    try
        [L_no, U_no] = lu_no_pivot(A);
        growth_no = max(abs(U_no(:))) / max(abs(A(:)));
    catch
        growth_no = inf;
    end
    
    % 选主元
    [L_yes, U_yes, P] = lu(A);
    growth_yes = max(abs(U_yes(:))) / max(abs(A(:)));
    
    fprintf('%-12s    %4d    %8.2e    %14.2e    %14.2e\n', ...
            matrix_names{i}, n, cond_A, growth_no, growth_yes);
end

%% 3. 部分选主元 vs 完全选主元
fprintf('\n3. 部分选主元 vs 完全选主元\n');
fprintf('目标：比较不同选主元策略的效果\n\n');

% 构造测试矩阵
rng(42);
A3 = randn(8, 8);
A3(1, :) = 1e-10 * A3(1, :);  % 第一行很小
A3(:, 1) = 1e-10 * A3(:, 1);  % 第一列很小
A3(1, 1) = 1;                  % 但(1,1)元素不小

fprintf('测试矩阵特点: 第一行和第一列很小，但A(1,1)=1\n');
fprintf('矩阵条件数: %.2e\n', cond(A3));

% 部分选主元（MATLAB默认）
[L3_partial, U3_partial, P3_partial] = lu(A3);
error3_partial = norm(P3_partial*A3 - L3_partial*U3_partial, 'fro');
growth3_partial = max(abs(U3_partial(:))) / max(abs(A3(:)));

fprintf('\n部分选主元结果:\n');
fprintf('重构误差: %.2e\n', error3_partial);
fprintf('增长因子: %.2e\n', growth3_partial);

% 模拟完全选主元的效果（简化实现）
[L3_complete, U3_complete, P3_complete, Q3_complete] = lu_complete_pivot(A3);
error3_complete = norm(P3_complete*A3*Q3_complete - L3_complete*U3_complete, 'fro');
growth3_complete = max(abs(U3_complete(:))) / max(abs(A3(:)));

fprintf('\n完全选主元结果:\n');
fprintf('重构误差: %.2e\n', error3_complete);
fprintf('增长因子: %.2e\n', growth3_complete);

%% 4. 迭代改进
fprintf('\n4. 迭代改进\n');
fprintf('目标：使用迭代改进提高解的精度\n\n');

% 构造病态线性系统
A4 = hilb(10);
x4_true = ones(10, 1);
b4 = A4 * x4_true;

fprintf('病态线性系统 Ax = b\n');
fprintf('矩阵: Hilbert矩阵 (10×10)\n');
fprintf('条件数: %.2e\n', cond(A4));
fprintf('真实解: x = [1, 1, ..., 1]''\n');

% 直接LU求解
[L4, U4, P4] = lu(A4);
x4_direct = U4 \ (L4 \ (P4 * b4));
error4_direct = norm(x4_true - x4_direct);

fprintf('\n直接LU求解:\n');
fprintf('解的误差: %.2e\n', error4_direct);

% 迭代改进
x4_improved = x4_direct;
fprintf('\n迭代改进过程:\n');
fprintf('迭代    残差范数      解的误差      改进倍数\n');
fprintf('----    --------      --------      --------\n');

for iter = 1:5
    % 计算残差
    r4 = b4 - A4 * x4_improved;
    residual_norm = norm(r4);
    
    % 求解修正方程
    delta_x = U4 \ (L4 \ (P4 * r4));
    
    % 更新解
    x4_improved = x4_improved + delta_x;
    
    % 计算误差
    error4_improved = norm(x4_true - x4_improved);
    improvement_factor = error4_direct / error4_improved;
    
    fprintf('%4d    %8.2e      %8.2e      %8.2f\n', ...
            iter, residual_norm, error4_improved, improvement_factor);
    
    if residual_norm < 1e-14
        break;
    end
end

%% 5. 扰动分析
fprintf('\n5. 扰动分析\n');
fprintf('目标：分析输入扰动对LU分解的影响\n\n');

% 基础矩阵
A5 = randn(8, 8);
[L5_orig, U5_orig, P5_orig] = lu(A5);

fprintf('扰动分析:\n');
fprintf('扰动水平    L的变化      U的变化      分解误差\n');
fprintf('--------    -------      -------      --------\n');

perturbation_levels = [1e-12, 1e-10, 1e-8, 1e-6, 1e-4];

for eps_level = perturbation_levels
    % 添加随机扰动
    rng(123);
    E5 = eps_level * randn(size(A5));
    A5_perturbed = A5 + E5;
    
    % LU分解
    [L5_pert, U5_pert, P5_pert] = lu(A5_perturbed);
    
    % 计算变化
    L_change = norm(L5_pert - L5_orig, 'fro') / norm(L5_orig, 'fro');
    U_change = norm(U5_pert - U5_orig, 'fro') / norm(U5_orig, 'fro');
    
    % 分解误差
    decomp_error = norm(P5_pert*A5_perturbed - L5_pert*U5_pert, 'fro') / ...
                   norm(A5_perturbed, 'fro');
    
    fprintf('%8.0e    %7.2e      %7.2e      %8.2e\n', ...
            eps_level, L_change, U_change, decomp_error);
end

%% 6. 条件数与稳定性
fprintf('\n6. 条件数与稳定性\n');
fprintf('目标：研究矩阵条件数对LU分解稳定性的影响\n\n');

% 生成不同条件数的矩阵
condition_numbers = [1e2, 1e4, 1e6, 1e8, 1e10, 1e12];

fprintf('条件数对稳定性的影响:\n');
fprintf('条件数      增长因子    前向误差    后向误差\n');
fprintf('--------    --------    --------    --------\n');

for target_cond = condition_numbers
    % 构造指定条件数的矩阵
    n = 10;
    [U_temp, ~, V_temp] = svd(randn(n));
    s = logspace(0, -log10(target_cond), n);
    A6 = U_temp * diag(s) * V_temp';
    
    actual_cond = cond(A6);
    
    % LU分解
    [L6, U6, P6] = lu(A6);
    
    % 增长因子
    growth_factor = max(abs(U6(:))) / max(abs(A6(:)));
    
    % 求解线性系统测试稳定性
    x6_true = randn(n, 1);
    b6 = A6 * x6_true;
    x6_computed = U6 \ (L6 \ (P6 * b6));
    
    % 前向误差
    forward_error = norm(x6_true - x6_computed) / norm(x6_true);
    
    % 后向误差
    residual = b6 - A6 * x6_computed;
    backward_error = norm(residual) / (norm(A6) * norm(x6_computed) + norm(b6));
    
    fprintf('%8.0e    %8.2e    %8.2e    %8.2e\n', ...
            actual_cond, growth_factor, forward_error, backward_error);
end

%% 7. 块LU分解的稳定性
fprintf('\n7. 块LU分解的稳定性\n');
fprintf('目标：分析块结构对稳定性的影响\n\n');

% 构造块矩阵
n = 12;
block_size = 4;

A11 = randn(block_size, block_size);
A12 = randn(block_size, block_size);
A21 = randn(block_size, block_size);
A22 = randn(block_size, block_size);

A7_block = [A11, A12; A21, A22];

fprintf('块矩阵结构: 2×2块，每块%d×%d\n', block_size, block_size);
fprintf('整体矩阵: %d×%d\n', n, n);
fprintf('条件数: %.2e\n', cond(A7_block));

% 标准LU分解
[L7_std, U7_std, P7_std] = lu(A7_block);
growth7_std = max(abs(U7_std(:))) / max(abs(A7_block(:)));

% 块LU分解（简化实现）
[L7_block, U7_block] = block_lu_decomposition(A11, A12, A21, A22);
A7_block_reconstructed = L7_block * U7_block;
error7_block = norm(A7_block - A7_block_reconstructed, 'fro') / norm(A7_block, 'fro');

fprintf('\n块LU分解结果:\n');
fprintf('标准LU增长因子: %.2e\n', growth7_std);
fprintf('块LU重构误差: %.2e\n', error7_block);

%% 8. 稳定性总结和建议
fprintf('\n8. 稳定性总结和建议\n');

fprintf('LU分解稳定性要点:\n\n');

fprintf('1. 选主元的必要性:\n');
fprintf('   • 部分选主元通常足够\n');
fprintf('   • 完全选主元在极端情况下更稳定\n');
fprintf('   • 增长因子是稳定性的重要指标\n\n');

fprintf('2. 条件数的影响:\n');
fprintf('   • 高条件数矩阵需要特别小心\n');
fprintf('   • 迭代改进可以提高解的精度\n');
fprintf('   • 后向误差通常比前向误差小\n\n');

fprintf('3. 实际应用建议:\n');
fprintf('   • 总是使用选主元的LU分解\n');
fprintf('   • 监控增长因子和条件数\n');
fprintf('   • 对病态问题考虑迭代改进\n');
fprintf('   • 使用双精度或更高精度计算\n\n');

fprintf('4. 数值稳定性判断标准:\n');
fprintf('   • 增长因子 < 10^6 (通常可接受)\n');
fprintf('   • 后向误差 ≈ 机器精度\n');
fprintf('   • 前向误差 ≈ 条件数 × 机器精度\n\n');

end

% 辅助函数
function [L, U] = lu_no_pivot(A)
% 不选主元的LU分解
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    if abs(U(k, k)) < eps
        error('主元为零，需要选主元');
    end
    
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U, P, Q] = lu_complete_pivot(A)
% 简化的完全选主元LU分解
n = size(A, 1);
P = eye(n);
Q = eye(n);
L = eye(n);
U = A;

for k = 1:n-1
    % 找到最大元素
    [max_val, max_idx] = max(abs(U(k:n, k:n)), [], 'all', 'linear');
    [i_rel, j_rel] = ind2sub([n-k+1, n-k+1], max_idx);
    i_abs = i_rel + k - 1;
    j_abs = j_rel + k - 1;
    
    % 行交换
    if i_abs ~= k
        U([k, i_abs], :) = U([i_abs, k], :);
        P([k, i_abs], :) = P([i_abs, k], :);
        if k > 1
            L([k, i_abs], 1:k-1) = L([i_abs, k], 1:k-1);
        end
    end
    
    % 列交换
    if j_abs ~= k
        U(:, [k, j_abs]) = U(:, [j_abs, k]);
        Q(:, [k, j_abs]) = Q(:, [j_abs, k]);
    end
    
    % 消元
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U] = block_lu_decomposition(A11, A12, A21, A22)
% 简化的2×2块LU分解
% A = [A11 A12; A21 A22] = [L11 0; L21 L22] * [U11 U12; 0 U22]

% 第一步：A11 = L11 * U11
[L11, U11] = lu(A11);

% 第二步：A12 = L11 * U12 => U12 = L11 \ A12
U12 = L11 \ A12;

% 第三步：A21 = L21 * U11 => L21 = A21 / U11
L21 = A21 / U11;

% 第四步：A22 - L21 * U12 = L22 * U22
S22 = A22 - L21 * U12;
[L22, U22] = lu(S22);

% 组装结果
n1 = size(A11, 1);
n2 = size(A22, 1);
L = [L11, zeros(n1, n2); L21, L22];
U = [U11, U12; zeros(n2, n1), U22];
end